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Abstract 

Various methods achieving importance sampling in ensembles of nonequilibrium trajectories enable one to 
estimate free energy differences and, by maximum-likelihood post-processing, to reconstruct free energy 
landscapes. Here, based on Bayes theorem, we propose a more direct method in which a posterior likelihood 
function is used both to construct the steered dynamics and to infer the contribution to equilibrium of all 
the sampled states. The method is implemented with two steering schedules. First, using non-autonomous 
steering, we calculate the migration barrier of the vacancy in Fe-a. Second, using an autonomous scheduling 
related to metadynamics and equivalent to temperature-accelerated molecular dynamics, we accurately 
reconstruct the two-dimensional free energy landscape of the 38-atom Lennard- Jones cluster as a function of 
an orientational bond-order parameter and energy, down to the solid-solid structural transition temperature 
of the cluster and without maximum-likelihood post-processing. 

Key words: Free-energy calculations, statistical thermodynamics, computer chemistry, molecular 
simulation 



1. Introduction 

One important application of molecular simulation is the estimation of the Landau free energy F of a 
given multi-particle system with respect to an order parameter £ 

F(Z) = -k B TlnP(£). (1) 

where T, ks and P{£) denote temperature, Boltzmann's constant and the probability to observe the system 
with value £ for the order parameter, respectively. Calculating Landau free energies thus amounts to 
measuring occurrence probabilities, a task that molecular simulation fails to achieve as soon as relevant 
portions of the phase space are rarely explored. So as to restore numerical ergodicity, many simulation 
techniques have been devised, based on umbrella sampling [1]. The generic idea of this technique consists 
in resorting to a judicious steering or restraining potential that enhances exploration of regions of phase 
space that would be poorly sampled otherwise. In its usual implementation, a series of umbrella sampling 
simulations [1| are first performed so as to cover the various regions of interest, and then the collected averages 
are combined using one of the various reweighing procedures @, HI, H| related to Bennett's acceptance ratio 
method [B| and based on likelihood maximization . 

In this context, Hummer and Szabo Q proposed to reconstruct the free energy profiles by applying the 
histogram reweighing procedure to nonequilibrium simulations Q instead of equilibrium simulations. To 
achieve this, they introduce an additional variable £ add and connect it to the relevant order parameter £ 
via the potential of umbrella sampling. Then, they mechanically steer the additional variable so as to push 
the particle system along the direction of the order parameter. They finally reconstruct the equilibrium 
properties by means of a two-step procedure. The first step provides the contribution to equilibrium at 
a given time-slice in trajectory space (after the collected nonequilibrium data have been reweighted using 
the probability ratios of the reverse-to- forward dynamics jjjl [lfj| within a path-average In a second 
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step, contributions arising from the entire range of times are finally combined using the weighted histogram 
analysis method 0, as in conventional umbrella sampling. 

We herein propose an estimator enabling one to retrieve equilibrium information from nonequilibrium 
trajectories. Like the aforementioned approaches based on Bennett's acceptance ratio method, our estimator 
resorts to reverse-to-forward probability ratios and retrieves information included in all time-slices. The 
estimator will have two advantages : (i) it does not involve any post-processing; (ii) it can be used with 
more general steering schedules than the one considered by Hummer and Szabo. We will illustrate these 
two points on the reconstruction of free energy landscapes. 

Concerning point (ii), we will consider Langevin dynamics in which steering arises from additional re- 
straining variables evolving stochastically and autonomously out of equilibrium into otherwise unexplored 
regions of phase space, enabling enhanced sampling along the steering directions. This way of proceeding 
can possibly be achieved by coupling the additional variables to high-temperature thermostats [12 . |l3L Il4 1 , 
as in multi-temperature sampling techniques 15, 16|], or by means of an adaptive biasing potential 17j . The 
former approach has been called temperature-accelerated molecular dynamics (TAMD) and the latter one 
metadynamics. 

The article is organized as follows. Section [5] establishes the general theoretical framework for the steered 
dynamics : the particle system and its additional variables are defined in subsection 12.11 the equations of 
the dynamics themselves are introduced in subsection 12. 2[ while the reverse-to-forward probability ratios 
associated with the dynamics, derived in subsection 12.31 are used to discuss the two steering schedules for 
the dynamics in subsection 12.41 In this framework, both the autonomous steering schedule of TAMD and 
the usual schedules that let a single steering variable evolve non-autonomously at constant speed Q appear 
as two particular limiting regimes. The reverse-to-forward probability ratios of subsection 12.31 are used to 
construct the two-state estimators 18, lj| 2(J 21, 22, 23 j reviewed in Section [3] and developed for calculating 
free energy differences with non-autonomous steering. Building on the approaches of Section [31 we propose 
in Section [4] an extended sampler and estimator enabling free energy reconstruction. The derivation that is 
given works both for autonomous and non- autonomous steering schedules. For completeness, we event ually 
show how the residence weight algorithm that is proposed relates to the waste-recycling algorithm [25j, |26( . 

Section [5] illustrates the performance of the proposed estimator with non- autonomous steering for a one- 
dimensional reconstruction problem. We compute the migration free energy of a vacancy in Iron (a-Fe). 
The second application given in Section [5] uses autonomous steering and reconstructs a two-dimensional 
free-energy landscape of the Lennard- Jones cluster with 38 atoms (LJ38). This benchmark system presents 
a rugged energy landscape, with two energy funnels separated by a high free energy barrier. It has been 
extensively studied using various methods, which eventually permits one to assess the relative numerical 
performance of our method with respect to existing methods. Concluding remarks are finally given in 
Section [71 



2. Extended Hamiltonian and steered dynamics 

2.1. Extended system 

Denote by r the particle position vector of dimension 3/ and by £ add an auxiliary vector of dimension 
J. The potential energy of the particle is E(r), while the steering potential is (l2l Il7l] 

where the order parameter £ = (£1, £j) °f dimension J is represented by a collective variable that is 
function of the particle positions. We denote the vector positions and vector momenta in the extended 
system by q and p, respectively. We have q = (r, £ add ) = (r 1; r 3 j, £f dd , £} dd ). The i-th components 
of these vectors are respectively denoted by qi and pi, with 1 < i < 31 + J. Let rm be the mass associated 
to the i-th component, and denote its momentum by pi — rriiqi where dots above coordinates designate 
time derivation. Denoting V(q) = E(r) + V(r, £ add ) the (total) potential energy of the extended system and 
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%(p,q) = Y^l=i J ^7" + V(q) its Hamiltonian, the normalized canonical probability density at temperature 
p- l =k B Tva 



pip, q) 



h 3I n 



Q /3.F-/3M(p,q) 



where the normalizing factor JF is the Helmholtz free energy of the extended system. It is related to the 
partition function logarithm 



T = -/TMn 



1 



h 3I Il 



dpdq 



Here, the infinitesimal volume with respect to coordinates reads 

3/+J 



dpdq = dpidqt 



(2) 



Canonical averages of any quantity A(r) defined with respect to particle positions can be taken in the 
extended ensemble (w denoting its phase space) as follows 



J A(r) exp [-PE(t)] dr 

J exp [-(3E(r)} dr 
J^(r)exp [-/3H(p,q)]dpdq 
L exp [-PH(p,q}\ dpdq 

A(r)p(p, q)dpdq 



(4) 



because contributions arising from the additional variables £ add can be inserted inside both integrals in 
Eq. [21 Hence, the additional variables and potentials do not affect the thermodynamic expectations of the 
particle system. 



2.2. Steered Langevin dynamics 

We consider that any coordinate qi with 1 < i < 31 + J is coupled to an independent thermal reservoir 
at temperature T. The traditional Langevin dynamics amount to propagating the system by solving the 
equations of motion (fi — —d qi H) 

q l =m~ 1 p i pi = fi- jiPi + b i (t)y / 2j. l k B Tmi (5) 

where bi represents a white noise of amplitude 1 and zero mean, while ji denotes the friction characterizing 
the coupling intensity with the i-th thermal bath. Here the amplitude of the fluctuations y/2yjk^Tml 
determines the temperature T. 

Let us assume that we have prepared the system in thermodynamic equilibrium at time t = [e.g. 
by propagating the Langevin dynamics ([5]) long enough and then setting the time to zero]. At t = 0, we 
switch on the external forces /J xt = ([ij — l)fj to act mechanically upon on the additional variables qj. The 
rescaling factors /ij are such that < [ij < 1 for 3/ < j < 31 + J. We also define /ij = 1 for 1 < i < 31 
by extension and we have /J xt + fj = fij fj . The frictional forces and the square of the fluctuations acting 
upon the additional variables are rescaled in the same way using the /ij's. The extended system is then 
propagated for a duration r using the steered Langevin dynamics below (using the convention on the indices, 
1 < i < 31 and 3/ < j < 31 + J) : 

qi = m^pi pi = fi - y t pi + b i (t)\/2~f t k B Tm l (6) 
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The rescaling of the dynamics amounts to unbalancing the interactions between particles and additional 
variables and to decreasing the coupling intensity with the j-th thermal bath while maintaining constant its 
temperature T. The effect of the external forces is to make the system depart from thermodynamic equilib- 
rium, by reducing the restraining forces acting upon the additional variables, which enhances exploration of 
phase space along the variables £ add . Additional commentaries about the dynamics (0 have been deffered 
to subsection 12.41 because they are based on the reverse-to- forward probability ratio derived in Section 12.31 
This ratio will indeed quantify the deviation of the dynamics with respect to equilibrium and will enable 
one to construct the estimators of Sections [3] and |H 



2.3. Reverse-to-forward probability ratio and discretization 

The probability to generate the dynamics along a path z = [q(t)]n< t < T given that the system is at 
(q(0),p(0)) at time t — is denoted Qf{z), while the probability to generate the same path using the 
reverse dynamics starting from system (q(T), p(r)) at time t — r down to t = is denoted Qr(z). Then, the 
reverse-to-forward probability ratio and the work W(z) done by the external forces fsf<j<3i+j u P on the 



extended system are related by the following expression (Eq. 2 in [37j |) 

P(p0),q0)) Qr{z) 



p(p(O),q(O))0F(«) 



= exp [-PW(z)] . (8) 



The identity above and the expression of the work will be derived explicitly for the discretized Langevin 
dynamics in this subsection. Before, we point out that identity ([5]) is similar to the more well-known 
identity involving the reverse-to- forward probability ratio due to Crooks [To} , except that a difference of free 
energy between a target system and a reference system appears in the latter form. No free energy difference 
appears here because the target and reference systems are the extended system itself with the extended 
Hamiltonian. The thermodynamic implications involving the two mentioned identities are compared in 
Ref. [37[. Besides, from a mathematical perspective, Eq. [8] can be interpreted as a generalized detailed 
balance equation involving the forward and backward Kolmogorov operators associated to our Langevin 



dynamics (see Eq. 4.43 in Ref. [24[). This interpretation allows both to define time reversibility rigorously 



and to extend the original derivations [l(J 37 1, which considered discrete-time Markov processes, to general 
continuous-time Langevin dynamics such as the one considered here. 

Since in practical applications we have to discretize the dynamics, we are authorized to expand the 
reverse and forward conditional probabilities, so as to include these quantities in the estimators directly. 
This is the approach that we follow in the sequel. Let At denote the discretization time step and Xn 
denote a state (q(i n ), p(t n )) at time t n = nAt. The discretized trajectory of a path z is characterized 
by the successive states (xo, ■ ■■,Xn, ■•■,Xn) obtained at times (to, ...,£„, ...,t./v) by propagating the Langevin 
dynamics forward star ting from a given state xo- This is achieved by updating the following discretization 
scheme [il [U [27|, [H, [29| from time t Q to time t N (1 < % < 31 + J) 

Pi,k+i/4 = Pi,ke~ 7tAt/2 + r)t,k+i/4 ( 9a ) 
Pi,k+i/2 = Pi,fc+i/4 + fi,kAt/2 (9b) 
Qi,k+i = qi.k + Pi, fe+i/2 At /mi (9c) 

Pt.k+3/4 = PiM+1/2 + fi,k+lAt/2 (9d) 

Pi,k+i = J3i,/c+3/4e~ 7,At/2 + (9e) 

where index k denotes time kAt while fi t k = ^ifi{kAt) and = (/i,; = 1 if i < 3/). Besides, the noises 
r ?j + /c+i/2±i/4 m Eqn. (|9a|) and (|9e[) are normal and have mean zero and variance o~i = (1 — e~^ iAt )mi/ (3. 
Updates (|9a|) and (|9e[) correspond to the momentum variations due to two consecutive Ornstcin-Uhlcnbeck 
processes of duration A* , These processes consist of propagating the momentum pi using 



Pi = ~liPi + bi(t)y/2miji/P 
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(10) 



from t = t n to t n + At/2 and from t = t n + At/2 to t n +i, where bi(t) is an uncorrelated white noise of unit 
amplitude. 

Because the scheme corresponds to a double Strang- Trotter decomposition [21[ with the position update 
in the center and then half momentum updates with respect to the force and the stochastic processes, it is 
symmetric and can thus be updated or downdated depending on whether the dynamics is considered to be 
forward or reverse. For the reverse dynamics, we must iterate 



Pi,k+3/± =Pi,k+ie ' + Vi, k +3/4 

Pi,k+l/2 = Pi,k+3/4 - fi.k+lAt/2 

q%,k = Qi,k+i - Pi, fc+1/2 At/rn, 

P»,fc+l/4 = Pi, fc+1/2 - /i,fcAt/2 

Pi,k = Pi,fc+i/4e _7iAt/2 + % k +i/4 



(11a) 

(lib) 
(11c) 

(lid) 

(lie) 



where the reverse noises r\ i k+i/2±i/4, nave the same variance. Besides, the time reversal of the Ornstein- 
Uhlenbeck process in (|lla[) and (|lle|) is the process itself. 

We then denote the probabilities of the discretized dynamics by P C ond(z|Xiv, N) and P C ond(z|xo: 0)- They 
will approximate the quantities Qr{z) and Qf{z) in (|24p. As a result of the discretization, the forward path 
probability can be factorized into the following product 



3/+JJV-1 



cond 



(z|Xo,0) 



n n $ -*(<fe+i/4)^(<fc+ 3 /4) 



i=l fc=o 
31+ J N-l 

n ii a 

i=l fc=o 



(2m,/?)- 1 



(^tfc+1/4) + (^tfc+3/4) 



(12) 
(13) 



where $ (Ti stands for the normal probability of variance <Tj = m, (1 - e~">" iA *)//3 and A ai denotes its nor- 
malizing factor. The normal laws are used to generate the stochastic noises vtk+1/4 anc ^ r ltk+3/4 °^ ^he 
i-th thermostat along trajectory z (0 < k < N). The conditional probability to generate z backward can be 
decomposed into a similar product of normal probabilities 



31+ J N-l 

n n ®°M,k+i/A)®°M k 

i=l k=0 
3I+.J N-l 

n n< ex p 

i=l fc=o 



i,k+l/il~* tT i yii,k+3/4 

(2m,/?)" 1 



1 - e -^ At 



(^^+1/4) + (^^+3/4)" 



(14) 
(15) 



Let Qi,k denote the temperature-scaled logarithm of the reverse-to-forward probability ratio associated 
with the Ornstein-Uhlenbeck processes for the i-th thermostat at step k. We have 



Qt 



(3 1 {in [$ CTi (^, fc+1/4 )<Mv +3/4 ) 



In 



i,fe+l/4 



] ^(vt, k +3/4> 



} 



(2m;) 



{ (^tfc+1/4) + (^tfc+3/4) (^^+1/4) + (^^+3/4) } 



1 / 2 2 1 2 2 \ 

\Pi,ft+l ~ Pi,k+3/4 + Pi,fe+l/4 ~ Pi,k J 



1 

2m, 



[Pi,fc+l -Pi,k] + 



At 



2 r 



r 2 _ f 2 

i,fc+l /i,fc 



:{li,k+l — <Zi,fe) ' \fi,k+l + fi,k)- 



(16) 
(17) 
(18) 
(19) 



The transformation from (|17[) to (|18[) involves expressing the noises as a function of the momenta after 
and before the Ornstein-Uhlenbeck processes and yields a form of detailed balance. In the transformation 
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from (|18[) to (|19p . the intermediate momenta Pi t k+i/4 an d Pi,fe+3/4 have been expressed as a function of the 
forces and positions at integer steps. 

The effective work done along the path from tg to t n = j^r defined by [2(| 



-r 1 



In 



P(X«) P cond(2:|Xri,«) 



conci 
3I+J n-1 



(20) 



i=l fc=0 



can be evaluated from the knowledge of the trajectory via the Qi,k's. The effective works will be used to 
retrieve equilibrium information in Section 01 Nevertheless, from a thermodynamical point of view, it is 
instructive to formulate the work W(z) in the continuum limit, achieved here when TV goes to infinity and 
with At = t/N. We thus define Q l (z) as 



N-l 



Q l {z) 



lim V" Qi 



N—>-\-oo 



(21) 



fe=0 



with At = t/N and all states x(^fc) inside the continuous path z when defining the limit. From the 
relation j20j 



JV-l 



Qi,k — 



k=0 

we deduce 



1 At 2 r - 

1 [2 2 1 | L ^ L 72 72 



1 



N-l 



~ o ^ (<?i,fe+l - <7i,fc) ' (/i,fc+l + /i,fe), 



(22) 



recalling that /i^/i = /i. Neglecting the constant Ito term arising from the integration of Pidpi, the quantity 
Qi(z) can be interpreted as the work done along the trajectory z by the force 

„ dp 4 



dt 



that is exerted by the i-th thermostat upon the i-th coordinate (-ji = /ij7i). This quantity thus represents 
the heat exchanged with the «-th thermostat. We recover an additional result in the continuum limit : the 
total heat exchanged with the thermostats during the forward dynamics, defined by 



3/+J 



Q(z) = &(*) 



(23) 



relates to the ratio of the reverse-to- forward conditional probability via the well-known expression 1^, II | 

Qr{z) 



Qf(z) 



exp [PQ{z)\ ■ 



(24) 



The heat also relates to the quantity W(z) defined in ([5]) via a conservation equation (l2"5j) . obtained by 
inserting ((24)) into (0 and then resorting to the relations p(p, q) = exp [(3 (F — H(p, q))] both at t = and 
t = t. We have 

Q{z) = H (p(r), q(r)) - H (p(0), q(0)) - W(z). (25) 



G 



Besides, resorting to fa = —d qi H in Eg . |2"21 and then summing yields an additional relation for the total heat 
exchanged with the thermostats 



= {H[p(T),q(r)]-H[p(0),q(0)]}- V (1 - fa)d qi Hd qi (t). (26) 



A 3I+J A 

i=l 



Note that the last summation runs from 3/ + 1 since fa = 1 for < i < 3/. Then, substracting [25] to 121)1 
enables one to identify W(z) explicitly 

3I+J , 9j (r) 

W(z) = / (l-MiJ^d^Ct). (27) 



j=3/+l 



?i(0) 



This quantity indeed corresponds to the work done by the external forces /J xt = (1 — p,j)d qj T-L upon the 
extended system, as stated in [301 ] . 

The limiting case consisting of setting fij = 1 for all j > 31 cancels the work in Eq. [57J [W(z) = 0]. This 
in turn implies a specific form of detailed balance p [x(t)] Qr{z) = p [x(0)] Qf(z) ensuring that the dynamics 
sample the equilibrium distribution p(x) =oc exp [—/3H (x)]- Aside from this limiting case, two particular 
schedules are possible for the steered dynamics depending on the /^--values in Eq. 1271 which we discuss 
below. 

2.4- Autonomous versus non- autonomous steering 

A first steering regime appears when the values of the scaling factors p,j 's are set to zero for all j > 31. The 
noise amplitude and the friction = pjjj vanish (see Eq. [7]). Any coordinate qj then evolves at a constant 
imposed velocity as in the schedule established by Hummer and Szabo The forces are conservative 
and time-dependent with respect to the real particles (once the additional variables have been eliminated 
by solving for them). The dynamics is said to be non-autonomous [3(j and we refer to this regime as 
non-autonomous steering. Non-autonomous dynamics with J = 1 guided by (time-dependent) conservative 
forces are well suited for computing free-energy profiles in one dimension or differences of free energy. Note 
that the fast switching schedule introduced by Jarzynski Q amounts to non-autonomous scheduling with a 
single external parameter X(t) = q3i+i(t) and p.31+1 = in Eg. 1271 Furthermore, the integral form in Eo. 1271 
with J = 1 corresponds to Jarzynski's definition of the work provided we consider the additional variable as 
a coupling parameter acting upon the Hamiltonian of the particle subsystem. 

The second steering regime consists of choosing < fj,j < 1 for j > 31. In this regime, the extended 
Langevin dynamics of Section |2"T21 is autonomous : the additional variables evolve stochastically by means 
of a force field {Mi/i}i<i<3/+j tnat i s time-independent and non-conservative [30| . i.e. that does not derive 
from a potential function except for particular conditions on the forces and the /itj's. We refer to this regime 
as autonomous steering. As will be shown in Section [S] autonomous steering is well adapted to the use of 
more than one additional variable. 

In the second regime, the dynamics of the additional variables may be given to a different thermodynamic 
interpretation. Indeed, the additional variables ([7]) also evolve according to the equation (fhj = p,J l mj, 
J > 3/) 



Qj 



h j ' /; - lAi + ; \ 'V'' 7 , 



where Tj — (J-J 1 T denotes the effective temperature of the thermostat that is actually coupled to qj and 
fhj denotes an effective mass for the additional variable. This dynamics is a particular implementation of 
the general dynamics given in [3lL Eq.3]. Dynamics coupled to thermostats at different temperatures reach 
a nonequilibrium steady-state with no well-defined temperature and satisfy a generalized detailed balance 
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equation 32|, |33| . Fluctuation theorems as well as reverse-to- forward probability ratios considered with 



respect to multi-temperature dynamics then relate to the heat transfers between the system and the various 
thermostats around the nonequilibrium steady state [32^ . The rescaling of the forces in the dynamics 
(subsection I2.2[) actually ensures that the reverse-to-forward path probability ratios relate to a transient 
mechanical work rigorously defined with respect to the equilibrium distribution of interest, as in the steering 
protocol of Hummer-Szabo. In particular, the work (|2"T)l depends on the potential energy of the extended 
Hamiltonian and not on its kinetic energy. 

Note that the stationary distribution reached by a multi-temperature dynamics exhibits a known an- 
alytical form 15, 1|1 12 1 when a separation of frequencies occurs between a slow variable qj subject to a 
thermostat at high temperature Tj and the remaining fast variables at normal temperature T. From this 
analytical form, the equilibrium probability profile P^? of qj at temperature T can be extracted from the 

established relation P!p(qj) oc [P st (qj)] Tj ^ T where P st (qj) is the stationary probability profile measured 
during a simulation. No separation of frequencies needs to be imposed in the approach of the present paper 
where the rescaling factor Tj/T = pj 1 acts upon the dynamics directly ©. 

In order to retrieve equilibrium information from transient nonequilibrium dynamics, Jarzynski derived 



its remarkable identity that involves the exponential average of the work (refer to [8|, |H, S S, HI HI for 
original and review papers on fluctuation theorems). We now briefly review the computational extensions 
that have been made to this approach for non-autonomous steering with a single additional variable. 



3. Two-state estimators for non-autonomous steering 



We assume here that trajectories are generated using a single steering variable A and with non-autonomous 
scheduling (//3/_|_i = 0). The phase space of non-autonomous paths is defined as follows 



n,. 



{z such that Vn, q3i+i{t n ) = X n } ■ 



A simple estimator associated to a biased sampler can give access to the ratio of normalizing constants 
related to the two thermodynamic states defined by Ao and A at [II Ell! El, Sil. This ratio can indeed 
be cast in the following form 



'3/+1 



'37+1 



Ao]p(x)dx 



/o na Pcond(z|xiv, N)p{xn)'Dz 

J Qaa PcondWXO, Q)p(xo)Vz 



J Qna [Pcond(s|XAT, N) P (xn)/Pb v (z)] Pb v (z)Vz 
Ta [Pcon< i (z\xo,0)p(xo)/PB ip (z)}P£(z)-Dz 



(28) 
(29) 



The first transformation f|28[) merely exploits the normalization of conditional probabilities with respect 
to path space f2 n a- The second transformation (|29[) formally inserts the biased probability distribution 
P v with respect to which sampling is performed (note that Jarzynski's identity is recovered by replacing 
P%(z) with P C ond(-g|yn 1 O)pn(Yn) where po(x) is the equilibrium density conditioned on (731+1 = Ao). In 
applications [H, HH [H, |39[ of identity (|29l) . the biasing potential tp is a function of the work function 
W(z) — —(3~ 1 \n{[P(z\xN, N)p(xn)] / [P(z\xq,0}p(xo)]}- In practice, the residence weight algorithm [l8| 
was observed to achieve good performance [19j, [23j (see Appendix O . 

The purpose of the paper is to extend the residence weight algorithm so that its sampler and associated 
estimator can handle the multiple thermodynamic states that can be defined owing to the extended system, 
irrespective of whether the scheduling of the steered dynamics is non-autonomous or autonomous. 



4. Multi-state estimator 

Residence weight algorithms can be formulated from two opposite points of view [23}. Here, we first 
build both the sampler and the estimator of the algorithm upon Bayes theorem by adopting the viewpoint 
of statistical inference. Then, we reinterpret the estimator as a conditional expectation (second viewpoint) 
in order to show the connection with the waste-recycling algorithm. 
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4-1. Posterior likelihood viewpoint 

Here, a marginal probability will be the importance function with respect to which path sampling is 
achieved, while a posterior likelihood function will be used on-the-fly to infer the equilibrium contribution 
of each generated state within the estimator. The marginal probability is defined in the path ensemble as 
the a priori probability of witnessing a path z under all possible hypotheses, i.e. as the sum of the product 
of all probabilities of hypotheses P v and corresponding conditional probabilities P CO nd: 

N 

P mW = E P cond WXn, ») ' P V (Xn, «). (30) 

71 = 

An hypothesis (xn, n ) is the knowledge of a state belonging to the path and of its index. The conditional 
probabilities in Eq. [3D] are given by 41, 27, 42| 



31+ J N-l 



Pcond(«|Xn,») = J] LI ^(^+l/4) $ -.(<fc+3/4) II (%k+l/d^i foi.fc+S/J- ( 31 ) 



-1 k—n k—n — 1 



where the normal distribution <$> ai are detailed in subsection 12.31 Here, P C ond(-z|Xri! n ) is the probability 
to generate the states Xn+i,Xn+2, ■■■Xn starting from any Xn by updating Eqn. I9all9el and then to gener- 
ate the states Xn-iiXn-2, •■■Xo starting from \ n by downdating Eqn. IllaTlllel The two particular cases 
Pcond(zl Yoi ) ~ Pf(z) and P C ond(z\xN, N) ~ Pr(z) were considered previously in biased path sampling 
schemes [20l |2~H I22I l39l . [23| to denote the probability to generate the forward and reverse trajectories, 
respectively. 

Additionally, the prior probability of hypothesis (x, n) in Eg. 1301 is 



P*(x,n) 



I P v (x)' lX ( n ) f° r non-autonomous scheduling, 
I P v (x) ivTT ^ or au t° nomou s scheduling. 



For non- autonomous scheduling, h x denotes the prior probability of index n and is such that h x (n) = 1 
if x pertains to the sliced phase space u) n = {(p, q) | <?3/+i = X n } that corresponds to index n, otherwise 
h x (n) — 0. Whatever x, we also assume that h x (n)h x (m) — for n 7^ m and ^2 71 =q h x (n) = 1 : the steering 
amplitude captures once all important regions of phase space. With autonomous scheduling, h x reduces to 
N , 1 whatever Xi involving the independence of the slice index from X- Besides, p v denotes a biased prior 
distribution of states 

P V (X) = exp - V { X ) - mix)] (32) 

where the biasing potential x ~^ fix) is here a state function (rather than a work function as in previous 
implementations) and the normalizing constant T v is the (^-dependent free energy. Irrespective of the 
scheduling, we have the useful equality 

JV . 

P V (X) = E / 5( X -Xn)P v (Xn,n)dxn- (33) 

where to denotes the unrestricted phase space (we previously assumed uj = Li^QU> n for non-autonomous 
scheduling). This equality will enable us to express the biased density p v as a path integral of marginal 
and posterior probabilities, irrespective of whether the scheduling is autonomous or non- autonomous. For 
this purpose, we introduce fi(xnj n) to denote the subspace of all paths going through Xn (at slice index n) 
and exploit the property that the sum over the path probabilities conditioned on Xn is normalized to one in 
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p v (x) 



N 

E 

n=0 
N 



1 concl 

{z\Xn,n)Vz 



Xn Jicond 

n JO 



E ~ Xn)Pscl{Xn, n\z) 



.n=0 



PUz)Vz 



(34) 
(35) 



where f a , n , Vz dxn simplifies to J Q Vz in Eq. [Ml with integration running over the space of either 

autonomous or non-autonomous paths (51 = 51 a or f2 na ). After permuting summation and integration in 
Eq. [MJ we introduced in Eq. [35]the posterior likelihood P sc i(x„,n|z) by resorting to Bayes relation 



Psd(Xn,n|z) 



A coiid 



(36) 



The posterior probabilities can be evaluated and simulated like the conditional probabilities. Indeed, plug- 
ging the various reverse-to-forward probability ratios given by 



Pcond(2:|X0,0)P y (X0,0) 
x concl 

(z\xn,n)Pv(xn,n) 



exp [<p(xo) ~ <p(Xn) ~ PW n 



into (|36|). yields the evaluable ratio 

Psei(Xn,n\z) 



exp [tp - (fin- PW n ] 

Ef=o ex P [<Po ~<Pk- pw k ] 



(37) 



(38) 



where W n is given in (1201) and ipk stands for (/?(Xfc) for simplifying. 

The aforementioned feature of the posterior and conditional probabilities makes it possible to generate 
a path distribution according to the marginal probability P^j(z). To explain how this can be done, let us 
consider a Monte Carlo move from Xn to x' n anc ^ whose associated transition probability Ptrans(XnlXn) obeys 
a detailed balance with respect to P v given by ([3"9"]) 



s(XnlXn)P ¥ '(Xn,n-) = Ptran S (Xn|Xn) P ^(Xn> n ) ■ 



(39) 



Then, considering a path z' containing (x' n , ri) and plugging Bayes relation (j3"6")l into ([59^1 for paths z and z' 
implies the detailed balance condition 

Pcond(z'|X™> «)Ptrans(XnlX™) P scl(X n , n\z)P^(z) = P C ond(z|Xn, ™)Ptrans (Xn |Xn)Psel(Xn, "I Z ')^M ( Z> ) ( 40 ) 
in which P cond (z'\x' n ,n)P tlli ns(x'„\Xn)Ps C \(Xn,n\z) and PcondWXn, »l)Ptran S (Xn|Xn)Psel(Xn' n \ z ') nave to be 

read as the probabilities to transit from path z to z' and from z' to z, respectively. The path distribution 
generated by any sampler satisfying the detailed balance condition 001 is ensured to converge toward the 
probability distribution Pj^. 

Still, the canonical average (HJ of quantity A(r) must be extended in order to be evaluable from a sample 
of paths distributed according to P^. To achieve this task, we first write the canonical average with respect 
to the biased probability measure. From the relation 

p{x) 



p v (x) 



exp \p (J--^) + V ?(x)], 
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we obtain 



(A) 



P V (X) 



p(x) 



p v (x)dx 



p v (x) 
■"°V(x)dx 



(41) 



Then, the path-integral expression of p v (Eq. l35j) is inserted into the biased average (Eq. |4"T|) and the Dirac 
functions are evaluated when integrating x over lo, which yields 



(A) 



N 

£ 

n=0 



A(r n )e^P scl (xn,n\z) 









r N 






P&{z)T>z 


1 


7 

Jn 


e^PsMn, n\z) 

- n=0 


P^Vz 


.(42) 



We are now in the position of evaluating the canonical ensemble average (|4"2"j) by following the traditional 



recipe 



i) we construct a Markov Chain distributed according to the probability density and ii) 



resort to the appropriate estimator to correct for the bias introduced by the importance function P^. 
(i) Repeat M times the following steps : 

(a) move to the shooting index n and state Xn of current path z which have been selected both with 
posterior likelihood P Be \(xn> n \ z ) m the previous steps (d) and (e); 

(b) generate the new shooting state x' n from Xn from probability Ptrans(Xralx™)> 

(c) set Co to 1, initialize the (provisional) next shooting index n prov to n and store x'n m t° Xprov! 

(d) shoot a new trajectory from state x'n an d index n (perform N — n updates from {x' n , n) followed 
by n downdates from (x' n ,n) again); concomitantly, at each Langevin iteration k, compute Cu = 



C h -i +exp 



m' n 



CO 



where index h = a 1 (k) runs from 1 to N after re indexing 

using a; with probability 1 — Ch-i/Ch, change n prov to a(h) — k and store x' a (h) m t° Xprovj 
otherwise leave n prov and Xprov unchanged; 

(e) set the next shooting index n' to 

Tlp r ov S-Ild Store Xprov 

into x'n' to denote the selected state of the 

new completed path z'; 

(f) go to (a) until the chain has been completed. 



(ii) Evaluate the estimator 



A 



t cxp(-flVf„| m ) 



M 



E"=qC*p(-'*'"I™-/ ?W "I™) 
E" = o°xp(-^„ h ) 

i:" = o°M-vn\ m -pw nlm ) 



(43) 



where {zi, z m , z^r} denotes the paths of the Markov chain constructed using the sampler and 



the simplified notations A r 



stand for A(r n ), ip(r n ) of path z m . W n i m represents the work 



done upon the extended system along the trajectory z m between xo and Xn (he. from times to to t r , 
via the mechanical coupling. 



The reindexing function is 
a(h) - 




if h < N - n, 
if h > N - n, 



k — n if k > n, 
N - k if k < n. 



Some details of the algorithm above such as the move from Xn to x' n m s t e P (i) _ (b) depends on the specific 
implementation: we choose Ptrans(XnlXn) = ^(x« ~ Xn) in Sections [5] implying that Xn is left unchanged; in 
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Section[5J x' n IS constructed from Xn by drawing new momenta p' n from the Maxwell-Boltzmann distribution. 
Besides, the estimator is obtained by plugging (1551) into the usual Metropolis estimator 

J Em=l En=0 ^»| m e y " |m P S cl(Xn| m , »|w|z m ) 

-Am = —j ]v ' i 

Em=i e v "i™P sc i(x„| m , n|m|z m ) 

related to ensemble average (jH]). Additionally, the shooting move [4l[ of step (i-d) generates the new path 
z' with probability P C ond(z'|Xn: n ) given in Eq. [31]by construction. The next shooting state yd , and next 
shooting index n' obtained from Xprov and n pmv are eventually selected with the compound probability 
(0 < h' < N) 



max(0,V-l) r _ _ 



n 



Ch Ch 



1 c " 




Ch' — Ch 



h=0 

where C_i = and ft,' = This compound probability is equal to (n' — a(h!)) 

Ch>-C h >-i_ exp[-<p' n , -W n >] 



Cn 



C » Ef= exp[-^-« 



As a result, the algorithm satishes the detailed balance equation P0|) . Note that the decomposition of the 
selecting procedure in (c-e) of (i) avoids storing all the configurations when a new path is constructed. 

4-2. Conditional expectation viewpoint 

In the residence weight algorithm, the shooting index n' related to path z' subsequent to z is constructed 
concomitantly with z' as outlined by steps (i)-(c) to (i)-(e) of subsection 14.11 As a result, the algorithm 
satisfies the following detailed balance condition 

P sel ( X ' n „n'\z')P cond (z'\x' n ,n)P^(x' n ,n) = P scl (x;,nk')Pcond(^|x^,«')P v (x^,«'), (45) 

for moves between (x'm n ) an d and also obeys detailed balance equation (1391) for moves between 

(Xn, n) and (x'n) n )- The algorithm thus leaves invariant the prior probability density P v : the Markov chain 
of states (xn , ri) constructed with the RW algorithm is thus distributed according to the probability density 

It follows from this theoretical description that the residence weight algorithm is a particular implemen- 
tation of the waste recycling algorithm introduced by Frenkel 2^, 23 1. Indeed, assuming there is no biasing 



potential in Eq. 1441 one can subsequently normalize the selection probabilities to one (M times) and write 

M N 

Am = jj ^2 ^2 A„| m P scl (x„| m ,n|m|z m ). (46) 

m— 1 n— 

which corresponds to the waste- recycling estimator given in Eq. 2.2 of Ref. [25[. Our symmetric selection 
procedure P so i corresponds to a Barker acceptance rule ([25, Eq. 2.1]) that considers states linked by the 
trajectories, while the wasted information included in the rejected Monte Carlo moves is recycled in the 
estimator of Eq. |4"61 

Interestingly, Delmas and Jourdain [44j showed that the estimator can be interpreted as the conditional 
expectation of A with respect to P se i and that it behaves normally asymptotically. Additionally, these 
authors proved that the statistical variance of the estimator is smaller than the one of the Metropolis- 
Hasting estimator when the acceptance rule is symmetric [44[ , as in the present situation. The first property 
implies that the conditional estimator is unbiased whatever the sample size M, at variance with maximum- 
likelihood estimators that are only unbiased in the limit of large samples |4[. The last property involving 
the variance reduction justifies the present strategy : when possible, one should systematically include the 
information contained in the states of the paths in statistical path-averages. 
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We now turn to the applications of our method. In Section [5j we implement non-autonomous steering 
along a one-dimensional order parameter with no biasing potential (ip = 0) and will resort to Eq. 05] to 
reconstruct the free energy profile. Then, in Section El the method will be tested on a more difficult 
benchmark model involving the direct reconstruction of a two-dimensional free-energy landscape. This task 
will be achieved by resorting to autonomous steering with respect to two additional collective variables. In 
addition, the biasing potential if will be constructed iteratively. 



5. ID free-energy reconstruction using non-autonomous steering 

The one-dimensional reconstruction problem involves calculating the migration free-energy of the vacancy 
in the a-Fe system. Atomic interactions of the model system are described by the (embedded atom method) 
potentials developed by Ackland and coworkers [45| and are computed using the minimum image convention. 
The crystal structure is body-centered cubic and the initial unrelaxed cell contains 1023 atoms displayed on 
1024 lattice sites : the vacant site (vacancy) is at a corner of the cell. Let ri = (91,92,93) denote a nearest 
neigbour atom of the vacancy along a [111] direction and define the one-dimensional collective variable 
dist(rx) as the distance between ri and the system's center of single additional variable 937+1 

is associated to dist(ri) via the potential energy §(937+1 — dist(ri)) 2 . The time-step is At = 4 ■ 10 _15 s and 
the friction of the Langevin dynamics is 7$ = 2.5 • 10 12 s _1 whatever i. 

Furthermore, protective spheres have been added upon the 7 nearest atoms of atom ri and upon the 7 
nearest atoms of the vacancy (ri being obviously unprotected). Each sphere is centred on the corresponding 
site of the underlying rigid lattice and is of radius a/2 (where a = 2.4728 • 10 _10 m is the nearest neigbour 
distance). Displacements moving a neigbouring atom out of its protective sphere are discarded, which can be 
done because the dynamics have been metropolized (see Appendix [B]) . The exit frequency of neighbouring 
atoms remains negligible even at high temperature. This procedure prevents spontaneous vacancy atom 
exchanges that may occur at temperatures above 540A" by the nonequilibrium steering without altering the 
statistics. 

The reaction coordinate £(r) is the projection of vector ri — R m along [111] direction. Measuring the 
quantity of interest, P(£), via a histogram amounts to monitoring the occupation probability of atom ri 
along [111] direction. We implement non-autonomous steering with its additional variable (£ add = 937+1) 
evolving at constant velocity according to 

937+l(*n) = 937+l(0) + - [q 3 I+l{t N ) - 937+l(0)] . (47) 
T 

The values 937+1(0) = —j^j and 937+1 (r) = in the steering schedule have been chosen such that the atom 
ri performs a single jump into the vacant site. The phase space is thus restricted to 2 possible vacancy sites. 
Time is given by t n = nAt. As noticed in related studies 4^, 47 1 implementing Jarzynski's work identity, 



we found it advantageous to use few long nonequilibrium trajectories rather than many short ones. Hence, 
for each temperature, we have generated M = 100 hybrid trajectories with N = 10 5 time-steps (implying 
a total of 10 7 force evaluations per simulation). From the 100 trajectories, we calculated the histograms 
Pioo(C) with £ spanning the interval [— A, ^] in 121 bins. Temperature ranges from 20K to 1000K. 

Figure [T] displays various outputs of the simulation carried out at T — 540K. Panel (a) displays the 
reaction coordinates £ of the states successively selected by P se i (Eq. I57]) . The high crossing probability 
discussed in the figure caption is related to the small value of the work performed on the system once the 
system has transited over the barrier at the present low-velocity steering. To illustrate the unavoidable lag 
effect caused by steering, let us consider the cologarithm of a Pi(£) estimate (i.e. of an estimate obtained 
from a single trajectory). The variation of — InPi(£) with respect to £ is related to the work done along 
the trajectory. We have represented in panel (b) the cologarithm of several Pi(£) estimates for forward or 
backward trajectories. The asymmetry of the colog-probability profiles is controlled by the steering direction 
and the amount of dissipation. The asymmetry is removed because the estimator combines trajectories 
generated forward and backward along £, i.e. starting from both free energy minima. The cologarithm of 
the Pioo(C) histogram, represented in panel (b) by the thick symmetric curve, illustrates this feature. A 
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Figure 1: Various simulation outputs expressed as a function of the reaction coordinate £ in meters. Panel (a) displays 
the trajectory number versus the § value of the selected states (x); panel (b) represents the colog of ten Pl(£) estimates 
measured from single paths generated either forward (red thin curve) or backward (green dashed curve); Panel (c) represents 
the divergence of the estimates. Since the biasing potential ip is zero, the selected reaction coordinates are distributed according 
to the equilibrium distribution P(£) itself. This one exhibits two maxima located by the position of the vertical dashed lines. 
The dashed segments that join the X symbols represent pairs of successively selected states. The fraction of dashed segments 
crossing the symmetric free energy barrier is about 40% at the considered temperature. 
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Figure 2: Free energy F(£,T) as a function of {; and T, estimated from — fcTlnPioo for the corresponding temperature. 

similar compensation has been observed previously (4cj | with a bidirectional variant of the Hummer-Szabo 
method. In panel (c), we have plotted the divergence defined by 



£>(£)=lnP(0-lnPi(0 



(48) 



where is estimated here using -Pioo(£) and the overbar denotes averaging the 100 available lnPi(£) 

values. Mathematically, the divergence -D(£) is a relative entropy betwen two distributions (49[. Thermo- 
dynamically, it is an excess entropy that is stored into the system and that would be irreversibly dissipated 
toward the thermostat if the system was allowed to relax back to equilibrium at constant £. This quantity 
gives information on the convergence of exponential averages (20L l49| . The smaller the divergence is and 
the more accurate the estimation is. We observe from FigfTJ that accuracy is smaller at the edges of the 
barrier, where the gradient of the steering potential is larger. Besides, -D(£) decreases again to a minimum 
around £ = 1.6 ■ 10~ 10 m that corresponds to the intermediate free energy minimum of panel (b). This trend 
suggests that the excess energy transiently stored in the steering potential is not entirely dissipated but is 
released to the extended system when £ reaches the intermediate energy minimum. 

The method was found to yield reproducible F(£, T)-estimates down to the temperature of 20K. Two 
series of simulations were carried out. From 20K to 250K, we used K = 10 • 10 4 J • mT 2 (IR1) and from 200K 
to 1000X we used K = 5 • 10 4 J • m~ 2 (IR2). Results are represented by the free energy landscape of Fig. [5J 
We observe that the intermediate free-energy minimum is more pronounced at the lower temperatures and 
completely disappears at temperatures higher than 700K . The migration free energies are deduced from the 
relative barrier heights along £-axis of Fig.[5J They are plotted as a function of temperature in Fig. [3] together 
with the prediction of classical harmonic approximation (CHA). CHA calculations have been performed using 
the procedure described in Ref. [271 . l50j and considering one of the two symmetric energy minima and saddle 
configurations. As expected, Monte Carlo simulations and CHA calculations agree at low temperatures 
(T < 200K) where anharmonic effects are negligible, confirming the exactness of our simulation method. At 
temperatures higher than 200-ftT, we observe a substantial deviation between simulation and CHA, attesting 
to strong an anharmonicity. Note that the extent of anharmonicity is in quantitative agreement with the 
one previously reported in the literature 51| for the vacancy migration free energy. 

In our first application, a single collective variable was used and simulations were performed successively 
with varying the temperatures so as to complete the landscape. In the second application, we show how 
to achieve two-dimensional reconstruction directly by resorting to autonomous steering with two additional 
variables. 
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Figure 3: Migration free energy as a function of temperature. IR1 and IR2 refer to two values of k (see text). The dotted line 
are the results of classical harmonic approximation. 



6. 2D free-energy reconstruction using autonomous steering 

We consider the 38-atom Lennard-Jones cluster. LJ^s is computationally troublesome to study because 



its potential energy landscape has two main funnels (52t l53l l54j , whose respective lowest energy structures 



are the icosahedron and the octahedron displayed in Fig. 2J It undergoes a two-stage phase change with 
increasing temperature starting from the octahedral structure. A solid-solid transition temperature between 
the octahedral funnel and the icosahedral funnel occurs near T ss sa 0.12e//cs, melting follows near Tr s ps 
0.17e/fcs. LJ reduced units of length, energy and mass (cr = 1, e/fcg = 1, to = 1) will be used in the 
following. 



As first collective variable, we use the bond-orientational order parameter Q4 of Steinhardt et al. [55 j 
that is able to distinguish between the cubic structures favoured at low temperatures and the icosahedral 
isomers above T ss . The second collective variable is the potential energy E{r). The associated additional 
variables, 53/+1 and (737+2, act upon the particles via harmonic potentials whose stiffnesses are k,\ — 10 4 
and K2 = 2, respectively. Their respective masses are 77137+1 = 6400 and 77737+2 = 0.8. The respective 



(a) , , (b) 



> 




Figure 4: The two lowest energy structures of the 38-atom cluster: (a) truncated octahedron with energy Eo = —173.9284 and 
(b) incomplete icosahedron, E\ = —173.2524. 
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Figure 5: Free energy F(Q 4 , E) as a function of Q 4 and E. Left panel is the actual measurement at T = 0.19, while the right 
panel represents the free energy reconstruction for temperature T = 0.05 as obtained by Legendre transform. 



coupling parameters are H31+1 = 0.9 and ^31+2 = 0.995, and the frictions are 7i = 5 ■ 10 3 [i < 31) and 
7j ; = 5 • 10~ 3 (j > 31). They have been chosen using the simple recipe that follows (j > 31) : (i) the fij's are 
tuned to enable the additional variables qj to oscillate with an amplitude large enough in the direction of 
the corresponding order parameter; (ii) the masses rrij are then tuned to set the velocity slow enough (but 
not too slow) and (iii) the coupling parameter jj is chosen small enough to ensure a smooth and regular 
evolution of the q/s. Procedures (ii) and (iii) prevent the dynamics from producing entropy, i.e. from 
dissipating the work done on the system into heat when the qj's evolve too fast. The values given above 
were found satisfactory and are certainly sub-optimal. Finding the optimal computational set-up is a non 
trivial task. 

A series of iterative simulations have been carried out at the temperature T = 0.19, using the procedure 
introduced by Coluzza and Frenkel [59j j in a similar context. Let P^(Q4, E) denote the histogram constructed 
by the l-th simulation. The biasing potential (p e+1 of the next simulation is then constructed using the 
iterative procedure 



1 - ^+1 



add 



ii+i 



(Q 4 ,E) = -\n(P i (Q i ,E)+p e n 



(49) 



The pf nin parameter determines the maximum value of the biasing potential and thus avoids possible sin- 
gularities arising from unexplored histogram bins. As the successive simulations explore larger portions of 
the phase space more and more accurately, the control parameter is decreased iteratively using the relation 
£min = 10 -9 ~ 4£ . Each simulation generates approximately M — 10 5 trajectories of N = 2.5 • 10 5 time-steps. 

Panel (a) of figure [5] represents the final Ft (E : Q 4 ) contour plot obtained at the temperature of the 
simulation To = 0.19. The color scale is such that the improbable regions of the landscape are displayed in 
black. The free energy at any temperature T± is related to the microcanonical density of states g{Q 4 , E) via 
the relation |56| 



F Tl (Q^E) 

k B Ti 



E 



r - In g(Q 4 ,E) + In J g(Q 4 ,E) exp 



E 



k B Ti 



dE. 



(50) 



Since this relation is also valid for the reference temperature To, it gives access to the difference of Landau 
free energies between any target temperature T\ and the reference temperature Tq of the simulation 



k B Ti 



E E-F To (Q 4 ,E) 



k B Ti 



k B T 



In / exp 



E-F To (Q 4 ,E) E 



k B T a 



k B Ti 



dE. 



(51) 



The integral over energy in (1511) acts as a normalizing factor and is the partition-function ratio involving 
the reference (To) and target (T\) systems. The free energy landscape at temperature T\ = 0.05 is fi nally 
displayed in panel (b) of Fig. [5] and reveals the low energy structures previously reported in Ref. 27 1. 
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To make a quantitative comparison between the present method [information retrieval with autonomous- 
steering, (IR)] and the three simulation methods used in Ref. [27| [nonequilibrium path-sampling (PS), 
parallel-tempering (PT) and Wang-Landau (WL)], we plot in Fig. [5] all the estimated free energy profiles, 
F(Q4, T) as a function of Q4, for temperatures T = 0.15, 0.12 and 0.05. 

For the present method and the Wang-Landau method [13] , we used the standard relation 



F{Q 4 ,T) = -k B T\n J exp 



'F T {Qi,E) 



dE (52) 



to obtain the free energy profile. In the reported path-sampling simulations, the temperature along the 
trajectories were slowly cooled down starting from T = 0.19, which, using the present terminology, amounts 
to non-autonomous steering with respect to temperature. However, in contrast with the present study, the 
estimator that was implemented was based on Crook's nonequilibrium average and could only exploit the 
information from the time-slice of the corresponding temperature. 

We observe that, at the temperature T = 0.15, the lowest that could be simulated correctly using 
standard umbrella sampling and histogram reweighting 1521 . our IR estimates for the free energies are in 
excellent agreement with the PS, PT and WL estimates [27] . However, at lower temperatures, we observe 
a disagreement for low free-energy structures in the range 0.02 < Q4 < 0.07, in particular for the one 
appearing around Q4 « 0.03, compared to estimates obtained using the three other methods (PS, PT and 
WL). This limitation of the method is con commitant to the slow convergence of the two-dimensional biasing 
potential <p a in this region of phase space. We did not perform another iteration as the total number of 



force evaluations was already 10 11 , i.e. equal to the one in the corresponding PS simulations |27[. 

Nevertheless, at temperature T — 0.05 [FigJBJ panel (c)], the estimates (IR) of F{Q^, T) revealed numer- 
ous low energy structures in the range 0.07 < Q4 < 0.12, which are explored with path-sampling, but missed 
with parallel tempering and Wang-Landau sampling. In particular, the low energy structure at Q4 — 0.12 is 
correctly predicted. Interestingly, the bassins of attraction of these minima are larger with IR than with PS. 
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Figure 7: Free energy F(Qg,T) as a function of Qg and T. The dotted line at T = 017 represents the phase transition between 
the liquid-like structure and the icosahedral structure (I). The one at T = 0.12 shows the transition between the icosahedral 
structure and the octahedral structure (O). The label D represents defected structures near Q6 = 0.25 (or Q4 = 0.08). 



This is due to the enhanced statistics made possible by the multi-state estimator that exploits information 
from all time-slices at any temperature, unlike the PS estimator. 

Finally, we plot the two-dimensional free-energy landcape F{Qq 1 T) as a function of Qg and T obtained 
from the simulation using a multi-state estimator in Fig. [71 The Qq order parameter is used because it 
better distinguishes the liquid phase L and the icosahedral phase /, as seen from this figure. Results are in 
qualitative agreement with previous simulation (PS) except for the defected structures D around Qe 0.24 
whose free energy is overestimated. 

To conclude this test study on LJ38 system, the reconstruction of a two-dimensional free-energy was 
achieved owing to the iterative construction of the two-dimensional biasing potential <j> with autonomous 
steering. In term of numerical efficiency, the approach was found advantageous because its estimator re- 
trieves the information contained in all the time-slices of the generated trajectories, unlike the nonequilibrium 
average previously implemented with non-autonomous steering [27]. Further comparing the numerical effi- 
ciency of the methods from Ref. [27| and the present one is difficult, as estimators and steering schedules of 
distinct types were used. The non-autonomous steering schedule of Ref. 27j should be tested with the more 
efficient multi-state estimator proposed in this study. In addition, given that the reported Wang-Landau 
sampling simulations achieved greater performance in the range 0.04 < Q4 < 0.08 but smaller performance 
in the range 0.08 < Q4 < 0.14, it may be worth implementing the multi-state estimator described in the 
present paper in combination with adaptive sampling methods such as Wang-Landau sampling 57, 58[ or 
metadynamics 17J so as to check whether the construction of the biasing potential would be facilitated. 



7. Concluding remarks 

In this article, we developed a unifying framework and a simple algorithm for retrieving the equilibrium 
information contained in all the time-slices from a sample of nonequilibrium trajectories. The algorithm, 
which shares several features with maximum-likelihood methods for nonequilibrium dynamics, is built upon 
Bayes theorem : a sampler generates a Markov chain of linked trajectories distributed according to a marginal 
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probability, while an estimator operates over the Markov chain so as to infer the contribution to equilibrium 
of the sampled data using a likelihood function. Contrary to maximum-likelihood estimators, the proposed 
estimator does not involve post-processing, can possibly be implemented with dynamics based on either 
non-autonomous or autonomous scheduling and is unbiased but not optimal. 

Concerning the overall efficiency of the method, we observe that, in agreement with theoretical pre- 
diction [40| | and maximum-likelihood simulations [48| . the most accurate estimations of free energies are 
obtained when trajectories can be initiated from the various regions of interest. This computational require- 
ment can be fulfilled by tuning an auxiliary biasing potential so as to flaten the prior probability distribution 
along the desired reaction coordinate. Using autonomous steering dynamics and resorting to a simple itera- 
tive procedure for constructing a two-dimensional biasing potential, the multi-state estimator could indeed 
reconstruct the free-energy landscape of the troublesome LJ38 system quite accurately. 

The presented simulations clearly outlined the advantages of steering autonomously rather than non- 
autonomously : the former strategy can be implemented in complex systems that require more than one 
steering variable and enables one to reconstruct multi-dimensional free energy landscapes directly. This 
feature in fact extends the possibilities of the Hummer-Szabo methodology [7J. 

Eventually, the multi-state estimator should be implemented with an adaptive sampler [g^, |6lj to check 
whether the combination of both techniques facilitates or not the construction of the biasing potential, 
compared to simulations resorting to one of the two techniques exclusively. This issue is to be considered in 
the wider perspective of waste-recycling [62|, 0| > which similarly advocates to retrieve all the information 
generated during the simulations within on-line statistical averages. A recent numerical investigation [r33| 
shows that a significant reduction of the statistical variances can be achieved in replica-exchange simulations 
that implement a multi-state estimator and a multi-proposal sampler 64] similar to the ones used in the 
present study. 
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A. Residence algorithms and optimized path-sampling 



In other implementations of the residence weight algorithm, the biasing function cf> is a work quantity 
rather than an auxiliary potential as in the present study. The equivalent form of our detailed balance 
equation (14TJ1) corresponds to the detailed balance condition [2J, Eq. 611 which itself formalizes the weighted 
balance condition given in earlier works (refer to [HI Eq. 5] and 0, Eq. 21]). The weighted balance 
equation was induced by analogy with the residence time algorithm. The latter algorithm [65l l66| and its 
extensions 67, 68, 6j| are used extensively in kinetic Monte Carlo simulations and also achieve importance 
sampling in ensembles of linked trajectories (more precisely, of kinetic pathways), owing to a similar selecting 
procedure. The selecting probability satisfies a detailed balance condition weighted by residence times (mean 
first passage times of exit). Eventually, the residence time algorithm also involves information retrieval [2^, 
[20! . I21I I23I ] . The kinetic pathways that can possibly be constructed from the master equation and that are 
eventually discarded by the algorithm [70| do contribute to the residence times. These analogies are more 
obvious for the residence weight algorithms used in Ref. 18, 1{| 2^, 21, 23l. which unlike the present case, 
generate and select trajectories pertaining to ramified paths called webs [711 ] . 

The relatively high numerical efficiency of the residence algorithms for estimating differences of free 
energies can be qualitatively explained by the study of Oberhofer and Dellago [40|. Assuming that the 
biasing potential is an adjustable functional depending on the work W(z), these authors derived the optimal 
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work dependent bias that leads to minimal statistical variance. The statistical variance was found to be 
minimal when the work-bias distribution contains typical forward and typical backward trajectories with 
similar weights, which is precisely an essential feature of the residence weight algorithms previously proposed. 
These algorithms generated paths using alternately the forward and reverse distribution in the extended 
ensemble of trajectories. 

In the present study, the biasing potential <p = cj> o £ is state dependent. The appropriate potential 
ensuring equipartition of trajectories would be such that </>(£) rj —/3F(£), which amounts to artificially 
flattening the probability density along the order parameter £. Optimizing the biasing potential ip requires 
knowledge of the quantities to compute, implying that the optimal bias can only be constructed iteratively, 
or adaptively, as outlined by Oberhofer and Dellago |40j . 



B. Metropolization of Langevin dynamics 

The unavoidable discretization errors are corrected in the ensemble average because the multi-state 
estimators uses the ratio of actual generating probabilities [2(J. Nevertheless, numerical efficiency depends 
on the choice of the time step At. Too small a time step decreases the sampling efficiency because states 
along the generated trajectories appear to be strongly correlated. Conversely, too large a time step produces 
numerical entropy (dissipated work pOj ) and one would observe that the selected states of the successive 
paths are also separated by small numbers of steps on average. Both situations results in increased statistical 
covariances in the path ensemble. Hence, in practice, the time step is tuned to achieve the best trade-off 
between decorrelation and entropy production. 

Unfortunately, for systems with many particles, discretization errors are important which imposes to 
choose a very small time step. This situation is encountered with our tabulated EAM potential of Iron. We 



have therefore metropolized the Langevin dynamics, i.e. we accept an iteration with probability [271 . I29j 

Pr = min(l, exp [-0 (H Xk (p k+1 , r fe+1 ) - H Xk ( Pfc , r fe ) - Q k )}) 

with the Metropolis rule. If the move is accepted, the new state is Xfe+i = (Pfc+i> r fc+i)> otherwise we set 
Xk+i = (— PfcjTjfc). The change of sign for momenta preserves the reversibility of the Markov chain and we 
take W{k — > k + 1) = H\ h+1 (xk+i) — H\ h {Xk+\)- A rejection rate of a few percent [I?} makes it possible 
to use a much larger time step, thus saving computational time. Note that when a rejection occurs with 
probability i_pa c c a ^ iteration k, this quantity must be included in the conditional probabilities P C ond(Xni n) 
whatever < n < N. Because the selecting probability P sc i considers ratios of conditional probabilities, the 
Metropolis rejections in the trajectories does not affect the work quantities W n and the algorithm given in 
subsection 14. II fsee fltA Appendix B.3] for the detailed proof). 

Finally, note that the two Ornstein-Uhlenbeck (OU) processes (j9aj) and (|9d|) of duration At/2 in the 
discretization scheme can possibly be merged into a single one of duration At. As a result, the discretization 
(Eqn. !9aH9dp , simplifies into a leap-frog scheme [2(| that generates a single noise per iteration 

Pi,k+i/2 = Pi,k-i/2 + Ui,k + fi,k) A< Qi,k+i = Qi,k + m~ 1 p i ^ +1 / 2 At. 

The quantity i^ k At describes the momentum variation during the OU process twice longer 
£ i>k At = (pi, k -i/2 + /i.fcAt/2) (e~^ At - 1) + Vhk = (f a -i/4 + kk+1/4) At/2 

where T]i t k is a normal noise of variance (1 - e- 2 ^ At )mi/p while £ hk „ 1/4 At/2 and £ hk+1/i At/2 denote the 
momentum variations during the OU processes in (|5af and (|9d|) . The leap-frog scheme could have been 
used in combination with the multi-state estimator since here the work does not depend on momenta at 
integer steps. 
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